%%   This script is the main file to replicate the results shown in Nucera, Sarno and Zinna, "Currency Risk Premiums Redux"  published in the Review of Financial Studies

%    In particular, the script produces: 
%    - Table 1, 2 and 3; 
%    - Fig. 1, Fig. 2 (table version along with the entries of Table A8 in the IA), Fig. 3 (along with FigA1 and FigA2 in the IA), Fig. 4 and Fig. 5 (table version) 
%   NOTES:
%   The script requries auxiliary functions written by others (see Read_me.txt for details). 
%   Full credit for these functions goes to Stefano Giglio, Craig Burnside, Anders Holtsberg and James P. LeSage.
    
%%
clear all
close all
clc

global xper xm2y nf nta T indc indcfm 

%% Add paths
addpath('FUNCTIONS');
set(0,'DefaultFigureWindowStyle','docked')
%% Test assets options 
w  = 20;   % weight on second pass [set equal to -1 to get back to PCA]
%% Set number of factors
pmin   =  1; % min number of factors 
pmax   =  6; % max number of factors       
%% Scaling factors
xper = 100;
xm2y = 12;
%% Giglio and Wu Options
indc     =  0; % 1 with constant (0 no constant)
%% Fama McBeth Options
indcfm   = 1;  % 1 with constant (0 no constant)
indcpe   = 1;  % if indcfm   = 1, 1 to include the constant in computation of the pricing error, 0 otherwise
%% ***************** LOAD DATA **********************
%% Load portfolio excess returns
load('DATA/R.mat')
%% Load G factors (candidate risk factors)
load('DATA/G.mat')
%% Load dates in matlab format
load('DATA/dateP.mat')
%% Load portfolios' labels
load('DATA/MatPrtLab.mat')
%% Load factors' labels
load('DATA/MatFctLab.mat')
%% String with fcts and pval
header  = ' ';
for i = 1:size(MatFctLab,1);
    header = strvcat(header, MatFctLab(i, :), '(pval)', '(se)');
end
header(1, :) = []; 
%% ***********************      MAIN CODE      ***********************
nf   = size(G, 2);
nta  = size(R, 2); 
%% Saving matrices [TO POPULATE]
MATmf = nan(nf, 5);
MATsrgc=nan(pmax,nf);
%% Loop over nf factors
for  j  = 1:nf
   
    strfctj = MatFctLab(j, :);
    
    %% annualize Factors and Test Asset Returns and remobe NaN
    Gj   = G(:, j);
    ifin = find(~isnan(sum([R Gj],2))); 
    r    = xm2y * R(ifin,:);
    g    = xm2y * Gj(ifin,:);    
    T      = size(g, 1);
    nwlags = round((4*(T/100)).^(2/9));
%% 2-pass Fama MacBeth (using ALL Test Assets)

iap  = 2;
max1 = nwlags;
max2 = nwlags;

tempindcfm=indcfm; 

[lamb, VO, VS, VG, Jmat, Rmat, ermod, b1] = gmm_fmb_lambdas(r, g, indcfm, max1, max2); 
[betatmp, pvaltmp]  = f_first_pass(r, g);
VX = VS; % Shanken       
lamb_g  = lamb;
lambt_g = lamb./sqrt(diag(VX));
lambse_g = sqrt(diag(VX));
lambp_g = f_pval2S(lambt_g, nta-(1+indcfm))';
jtest   = Jmat(2, 2);    

if indcpe ==0
mae     = abs(ermod(:,1)'-ermod(:,2)');
else
mae     = abs(ermod(:,1)'-ermod(:,end)');
end  

MATbeta_g(:, iap, j)   = betatmp(:, 2);
MATbetap_g(:, iap, j)  = pvaltmp(:, 2);
MATlam(j, iap)         = lamb_g(1+indcfm); 
MATlamp(j, iap)        = lambp_g(1+indcfm); 
MATlamse(j, iap)       = lambse_g(1+indcfm); 
MATrp_g(:, iap, j)     = ermod(:,2)';
MATmae(:, iap, j)      = mae(:); 
MATjtest(j, iap)       = jtest; 

%% 3-pass Giglio and Wu (2019)    
iap  = iap+1:iap+pmax;
r_gw =  xm2y*R;
g_gw =  xm2y*Gj;

res  = rppca_rfs(r_gw', g_gw', pmax, nwlags, indc, j, w); 

ifinG = find(~isnan(sum(Gj,2))); 
if (size(ifinG,1) == size(dateP,1)) %save vhat,bhat,gammatilde, pricing performance, pmax for the full time span  
 vhat=res.vhat; %save vhat for later use
 betahat=res.betahat;%save betahat for later use
 MATmSR=res.MaxSR;       
 MATrms=res.RMS;            
 MATsigmae=res.sigmae;       
 MATweights=res.weights;
 MATiSR=res.iSR;
 MATmaxnf=res.maxnf;
 MATeig = res.eigenvalues;
 Gammatilde=res.gammatilde;%save latent factor prices of risk for later use
 pe_tilde = res.alphahat; %save latent factors associated pricing errors.
 R2_lf    = res.R2F; %save latent factors R2.
 pmaxes   =res.pmax; % save maximum number of latent factors suggested by the test
end
 MATlam(j, iap)        = res.Gammahat(1+indc, :)'; 
 MATlamp(j, iap)       = res.Gammahatp(1+indc, :)';
 MATlamse(j, iap)      = res.Gammahatse(1+indc, :)';
 MATr2g(:,j)           = res.R2G;  
 MATwtest(j,iap)       = res.wstat';  
 MATwtest(j,1:2)       = NaN;
 MATwtestp(j,iap)      = res.wtestp';
 MATwtestp(j,1:2)      = NaN;
 %% Compute the (cleaned) g-factor Sharpe ratio 
for nv=1:pmax
gclean=res.etahat(:,1:nv)*vhat(:,1:nv)';
gclean=gclean';
MATsrgc(nv,j)=abs(sqrt(12)*nanmean(gclean)./nanstd(gclean));
end

end
%% store pricing performance of latent factor model
MATperfll_mv = [R2_lf'.*100  mean(abs(pe_tilde),1)'.*100];                  
%% Make Tables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')  
disp('*****Table 1')
disp(' ') 
infoT.fmt  = '%10.2f';
infoT.rnames = strvcat(strcat('$\gamma$= ',num2str(w)),strcat('F(',num2str(1),'-',num2str([1:pmax]'),')'));
infoT.cnames = strvcat('sig_e', 'RMS', 'R2(%)', 'MAE(%)', 'SR', 'b_mv_k', 'mu_f'); 
tabtmp = [100*MATsigmae'  100*MATrms' MATperfll_mv(:,1) MATperfll_mv(:,2)  sqrt(12)*MATmSR'   MATweights(end,:)' mean(100*vhat)'];
mprint(tabtmp,infoT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

LabGFctPval  = ' ';
LabGFctexp   = ' ';
for j = 1:nf
TAB_lam(j*5-4:j*5, :) = [[MATmf(j, 1) MATlam(j, :)];...
                         [MATmf(j, 2) MATlamse(j, :)];...
                         [MATmf(j, 3) MATwtestp(j, :)];...
                         [MATmf(j, 4) [NaN(1,3) MATsrgc(2:pmax,j)']];...
                         [MATmf(j, 5) [NaN(1,3) 100*MATr2g(2:pmax,j)']]];

LabGFctexp  =  strvcat(LabGFctexp  , MatFctLab(j, :));                     
LabGFctPval  = strvcat(LabGFctPval , MatFctLab(j, :), '(se)', '(pval_weak)', 'SR', 'R2');
end
LabGFctexp (1, :) = []; 
LabGFctPval (1, :) = []; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')  
disp('*****Table 2')
disp(' ') 

% loop over g factors to recover eta estimates, standard errors and pvals
mattab2 = [];
for  j  = 1:nf
FF   = vhat;
kk   = pmax;
GGj  = xm2y * G(2:end, j);
iFin = isfinite(GGj);
GGj  = GGj(iFin, 1); 
FFj  = FF(iFin, :);
nlag = round((4*(sum(iFin)/100))^(2/9));
YY   = GGj;  
XX   = [ones(size(YY)), FFj(:,1:kk)];
output  = nwest(YY, XX, nlag);
etatmp =  output.beta';
etase  = (output.beta./output.tstat)';
etap   = f_tstat2pval(output.tstat, sum(iFin)-kk);
    
% get R2 by sdf dimension
for isdf = 1:pmax
    
YY          = GGj;  
XX          = [ones(size(YY)), FFj(:,1:isdf)];  
output      = nwest(YY, XX, nlag);
etaisdf     =  output.beta';
etaseisdf   = (output.beta./output.tstat)';
etapisdf    = f_tstat2pval(output.tstat, sum(iFin)-isdf);
etar2       = output.rbar;

matetar2isdf(j, isdf, 1)  = output.rsqr; 

end

matetaols(j, :, 1) = etatmp(2:end); 
matetase(j, :, 1)  = etase(2:end); 
matetap(j, :, 1)   = etap(2:end); 
end

% to display etas and R2
nfoT.fmt  = '%10.2f';
infoT.cnames = strvcat(strcat('eta_F',num2str([1:pmax]')), strcat('R2_F',num2str([1:pmax]')));
infoT.rnames = strvcat(' ', MatFctLab(1:nf, :)); 
tabgfctexp = matetaols;
tabgfctr2  = xper*matetar2isdf;
mprint([tabgfctexp tabgfctr2], infoT);   

% to display pvals of eta estimates (to get stars)
nfoT.fmt  = '%10.2f';
infoT.cnames = strvcat(strcat('pval_F',num2str([1:pmax]')));
infoT.rnames = strvcat(' ', MatFctLab(1:nf, :)); 
tabgfctpval = matetap;
mprint([tabgfctpval], infoT); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' ')  
disp('*****Table 3')
disp(' ') 

TAB_lamP = TAB_lam;
infoT.fmt  = '%10.2f';
infoT.cnames = strvcat(strcat('FMB(', num2str(size(r, 2)),')'), strcat('F(1-',num2str([2:pmax]'),')'));
infoT.rnames = strvcat(' ', LabGFctPval); 
mprint([TAB_lamP(:, [3 5:end])], infoT)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Make Fig 1: Tradeoff
load('DATA/hmldata.mat')
NP = 5;
Fig_1

%% Make Input for Fig 2 and Table A8 (plain and normalized eigenvalues) 
Rnan=R; % drop gap filled observations
Rnan(389:411, 42:46)=nan;
[tabtmp, txteig, vare] = f_get_eigenvalues(vhat, betahat, Rnan, w, pmax);

%% Make Fig 3:  Portfolios risk exposures to latent factors
load('DATA/Fig3_input', 'TypTA', 'MatNP', 'LabPrtSel', 'LabCSSel', 'LabTSSel', 'TypCSSel', 'TypTSSel', 'pmax', 'optordpc')
[matr2prt] = f_fig3(betahat, vhat, Rnan, TypTA, MatNP, LabPrtSel, LabCSSel, LabTSSel, TypCSSel, TypTSSel, pmax, optordpc);

%% Make Fig 4: Observed vs Model-Implied Portfolio Excess Returns
Fig_4
%% Make input for Fig 5 and plot and example: Nontradable Factors’ Explained Variation by Latent Pricing Factors
fprintf(['\n'])
fprintf('Figure 5 (table version): Nontradable factors’ explained variation by latent pricing factors\n')
fprintf(['\n'])
tab = MATr2g(1:3,:)';
infoT.fmt  = '%10.2f';
infoT.cnames = strvcat(strcat('R2(%); F(1-',num2str([1:3]'),')'));
infoT.rnames = strvcat(' ', LabGFctexp); 
mprint(100*tab, infoT);

%example nfsel g factor (not ordered)
figure('Name', 'Fig5 (example)') 
nfsel =10; %must be <=nf
bar(100*[tab(1:nfsel, 1) tab(1:nfsel, 2)-tab(1:nfsel, 1) tab(1:nfsel, 3)-tab(1:nfsel, 2)],'stacked')
xticklabels(MatFctLab(1:nfsel, :))
ylabel('R2(%)')


